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Abstract: In this paper we present deflation and augmentation techniques that have been de- 
signed to accelerate the convergence of Krylov subspace methods for the solution of linear systems 
of equations. We review numerical approaches both for linear systems with a non-Hermitian coef- 
ficient matrix, mainly within the Arnoldi framework, and for Hermitian positive definite problems 
with the conjugate gradient method. 
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Techniques de deflation et d'augmentation dans les solveurs 

lineaires de Krylov 



Resume : Dans ce rapport nous presentons des techniques de deflation et d'augmentation 
qui ont ete developpees pour accelerer la convergence des methodes de Krylov pour la solution 
de systemes d'equations lineaires. Nous passons en revue des approches pour des systemes 
lineaires dont les matrices sont non-hermitiennes, principalement dans le contexte de la methode 
d'Arnoldi, et pour des matrices hermitiennes definies positives avec la methode du gradient 
conjugue. 

Mots-cles : Augmentation, Deflation, Methodes de Krylov, Systemes lineaires d'equations, 
Preconditionnement . 
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1 Introduction 

The solution of linear systems of the form Ax* = b plays a central role in many engineering and 
academic simulation codes. Among the most widely used solution techniques are the iterative 
schemes based on Krylov subspace methods [3J (23J [Ml HO]- Their main advantages are their 
ability to solve linear systems even if the matrix of the linear system is not explicitly available 
and their capability to be "easily" parallelizable on large computing platforms. In order to speed 
up the convergence of these solution techniques, Krylov subspace methods are almost always 
used in combination with preconditioning. That is, instead of solving directly Ax* = b, the 
linear system is transformed into an equivalent one, e.g., Mi Ax* = M\b, referred to as left 
preconditioned system, that is expected to be more amenable to a solution. The definition of an 
efficient preconditioner M\, that should be an good approximation of A -1 in some sense, is very 
much problem dependent and is consequently an extremely active research field. We can also 
consider other equivalent linear systems AM^t* = b with x* = M^t* (right preconditioner) or 
M\AM2t* = M\b with x* = M%t* (split preconditioner). We refer the reader to [8] for a detailed 
overview on preconditioning. 

There exist two complementary alternatives to speed up the convergence of the Krylov space, 
namely augmentation and deflation. Roughly speaking, in augmentation techniques, the search 
space in an enlarged Krylov space that is defined by a direct sum of two subspaces. This search 
space Si (of dimension £) has the following form 

S e = K m {A,b)®W (1) 

where /C m (j4, b) is a Krylov subspace of dimension m generated by the matrix A and the vector 
b and W (of dimension k) is called the augmentation space. A typical goal of augmentation is to 
add information about the problem into the global search space Se that is only slowly revealed 
in the Krylov subspace itself. 

Alternatively, deflation is based on the use of a projection operator P to decompose x* as 
x* = Px* + (I — P)x* . The general idea is to select P such that the solution of PAx* = Pb, 
referred to as the deflated linear system, is more easily amenable to a solution by a Krylov 
subspace method than the original linear system Ax* = b. The component (I — P)x* can then 
be computed by solving a linear system of small dimension. 

The purpose of this paper is to expose these two latter acceleration techniques that become 
increasingly popular. We refer the reader to [3D] for a recent excellent analysis of these methods 
together with detailed references and historical comments. Here augmentation and deflation are 
described in a framework where variable preconditioning can be used as it is nowadays customary 
when considering large scale linear systems [HU HH E3] ■ This paper is organized as follows. In 
Section [2] we introduce some background on Krylov subspace methods with emphasis on the 
minimum residual norm approach for systems with a non-Hermitian coefficient matrix and the 
conjugate gradient method for the solution of Hermitian positive definite problems. In Section [3J 
we describe the augmentation and deflation techniques and their possible combination in the 
case of systems with non-Hermitian matrices with references to concrete applications. Similar 
exposure is performed in Section [¥] for Hermitian positive definite linear systems. Finally some 
concluding remarks and prospectives are drawn in Section [5J 

2 Some background on Krylov subspace methods 

We briefly describe the basic properties of Krylov subspace methods for the solution of a linear 
system of equations of the form 

Ax* = b (2) 
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where the nonsingular n x n coefficient matrix A € C" x ™ is supposed to be either non-Hermitian 
or Hermitian positive definite and b a given vector in C". First we introduce the notation used 
throughout this paper. 

2.1 Notation 

We denote the range of a matrix A by Rangc(A) and its nullspace as Ker(A). We denote 
by |.| the Euclidean norm, Ik G C kxk the identity matrix of dimension k and Oi X j € C IX - J 
the zero rectangular matrix with i rows and j columns. T denotes the transpose operation, 
while H denotes the Hermitian transpose operation. Given a vector d £ C fe with components 
di, D = diag(di, . . . , d/j) is the diagonal matrix D 6 <C kxk such that Da = d{. Given Z m = 
[zi, ■ ■ ■ , z m ] € C nxm we denote its i-th column as z t € C™, (1 < i < m). The vector e m e E m 
denotes the m-th canonical basis vector of R m . We denote by k(A) the Euclidean condition 

number of A that is defined by n{A) = max where cr max (<r m in) is the largest (respectively 

^min 

smallest) singular value of A. For Hermitian positive definite matrices, the condition number 

reduces to k(A) = ^ max where A max (A m ; n ) is the largest (respectively smallest) eigenvalue of 

A. Finally, throughout the paper for the sake of readability the integer subscript I denotes the 
dimension of the search space. 

2.2 Basic properties of Krylov subspace methods 

The Krylov subspace methods seek for the solution of Equation ^ in a sequence of embedded 
spaces of increasing dimension K.i(A,b) = span(b 7 Ab, ■ ■ ■ This is motivated [44] by 

the fact that for £ large enough these spaces contain the solution of the linear system If 
we denote by m^x) the minimal polynomial associated with A, the Jordan decomposition of 
this polynomial writes mA(t) = — K) mi where (Ai, • • • , A s ) are the distinct eigenvalues 

of A and (mi, • ■ ■ ,m s ) their indices in the Jordan form. In a canonical form we also have 
mA^t) = YliLo a it l with m = Si=i m i an d a a = Y[i=i(^^i) mi 7^ since A is nonsingular. 
Consequently, A^ 1 = — a^ 1 o 1 a i+iA l that portrays x* = A~ x b as a vector of the Krylov 
space K, m -i(A, b). This indicates that, in exact arithmetic, Krylov methods must converge in at 
most m — 1 steps or less if the right-hand side does not have components in all the eigendirections. 
This observation also gives some ideas on ways to speed-up the convergence of these methods. 
As mentioned earlier preconditioning is a widely used approach that consists in transforming @ 
in an equivalent nonsingular system where the preconditioned matrix has less |58j or better 
clustered eigenvalues (see [5] and the references therein). 

The rest of the paper is dedicated to an overview of proposed techniques for augmentation 
and deflation both for non-Hermitian and Hermitian positive definite problems. 

2.3 Minimum residual Krylov subspace method 

In this section we focus on minimum residual norm subspace methods for the solution of linear 
systems with a non-Hermitian coefficient matrix. We refer the reader to [Ml [HO] f° r a general 
introduction to Krylov subspace methods and to [73] for a recent overview on Krylov subspace 
methods; see also [23 [21] for an advanced analysis related to minimum residual norm Krylov 
subspace methods. 

Augmented and deflated minimum residual norm Krylov subspace methods are usually char- 
acterized by a generalized Arnoldi relation introduced next. 
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Definition 1 Generalized Arnoldi relation. The minimum residual norm subspace methods in- 
vestigated in this paper satisfies the following relation: 

AZt = V e+1 H e (3) 

where Z e G C nxe , V e+1 G C ,IX ^ +1) such that V^Ve+i = h+i and H e G C {e+1 ^ xe . These methods 
compute an approximation of the solution of @) in a i-dimensional affine space xq + Zi yi where 
yt G C . In certain cases, Hi is an upper Hessenberg matrix. 

We next introduce a minimum residual norm subspace method proposed by 
Saad [57] since it is the basis for further developments related to augmented and deflated Krylov 
subspace methods of minimum residual norm type. This method named Flexible GMRES (FGM- 
RES) was primarily introduced to allow variable preconditioning. We denote by Mj the nonsin- 
gular matrix that represents the preconditioner at step j of the method. Algorithm [JJ depicts the 
FGMRES(^) method where the dimension of the approximation subspace is not allowed to be 
larger than a prescribed dimension noted I. Starting from an initial guess xq G C n , it is based 
on a generalized Arnoldi relation 

AZ t = Vt+iH t with V^Ve+i = h+i, (4) 

where Zi G C nxe , Vi+i G C nX ^ +1 ' and the upper Hessenberg matrix Hi G C^ +1 ^ xf are obtained 
from the Arnoldi procedure described in Algorithm [5J An approximate solution xi G C" is 
then found by minimizing the residual norm ||6 — A(x$ + Ziy)\\ over the space xq + range(Z^), 
the corresponding residual being ri = b — Axi G C" with ri G range(V^+i). With notation of 
Algorithm [T] the current approximation xi can be written as 

xi = x Q + Ziy*, (5) 

whereas the residual ri — b — Axi satisfies the Petrov-Galerkin orthogonality condition 

n _L A range(Z^). 

Hence, an optimality property similar to the one that defines GMRES is thus obtained [5H] . We 
note however that no general convergence results are available since the subspace of approximants 
Range(Zf ) is no longer a standard Krylov subspace. We refer the reader to [671 R)5] for the analysis 
of the breakdown in FGMRES. Furthermore, as it can be seen in Equation the update of the 
iterate xi requires to store the complete set of vectors Zi inducing a large memory footprint for 
large I. In order to alleviate this memory requirement, a restarting strategy must be implemented 
as shown in Algorithm [TJ The construction of a complete set of Zi is often name a cycle of the 
method and corresponds to one iteration of the loop in Algorithm [T] 

When the preconditioner is constant, FGMRES (£) reduces to right-preconditioned GMRES (£) 
whose convergence properties are discussed in [SSI Chapter 6] . 

2.4 Conjugate gradient method 

The conjugate gradient [IT] is the method of choice for Hermitian Positive Definite (HPD) linear 
systems. In a shortcut, it relies on an Arnoldi like relation (namely a Lanczos relation |51| ) 
similar to Equation Q with Zi = Vi (case of no preconditioning) and a Ritz-Galerkin condition 
ri = b — Axi J- Ki(A,ro). At each iteration xi = xq + Viy* is computed via the solution of 
the small linear system Hiy* = ||ro||(l, Oix(i-i)) T > where Hi = V e H AVi is the square leading 
part of Hi. Since A is Hermitian, Hi is also Hermitian. Furthermore, its structure is upper 
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Algorithm 1 Flexible GMRES(f) 

1: Initialization: Choose £ > 0, tol > 0, xq S C". Let tq — b — Axq, (3 = \\ro\\, c = [j3, 0i x e} T 
where c G C i+1 , V\ = r /f3. 
Loop 

2: Computation ofVe+i, Ze and Hi (see Algorithm^): Apply £ steps of the Arnoldi method 
with variable preconditioning (zj = MJ Vj,l < j < £) to obtain Ve+\ G C x ' £+1 ', 
Ze G <C nx( and the upper Hessenberg matrix He G C" +1 ' xt such that 

AZ e = Vt+i&i with V^Vi+i = h +1 . 

3: Minimum norm solution: Compute the minimum norm solution X£ G C" in the affine 
space xq + range (Z^); that is. xt = xq + Ziy* where y* = argmin ||c — H£y\\. 

4: Check the convergence criterion: If ||c — -ff^?/*||/||6|| < tol, exit 
5: Restarting: Set a; = xt, r = b - Ax , (3 = \\r \\, c = [/3, 0i x e] T , V\ = r //?. 
End of loop 



Algorithm 2 Arnoldi procedure: computation of Ve+i, Z( and He 
l: for j = l,i do 

2: Zj - = Mj\ 
3: S = Az 3 ' 

4: for i — l,j do 

5: hij = V^S 

6: S = S - hi,jVi 

7: end for 

8: h i+ ij = ||s||, = 8/hi+lj 

9: end for 

10: Define Zi = \z\, • ■ ■ , ze], V e+ i = [vi, ■ ■ ■ ,Vi+i], H e = {hij}i<i<i+i t i<j<i 



Hessenberg that, combined with the Hermitian property, implies that He is tridiagonal HPD. 
The first consequence of this structure of He is that the orthogonalization of Ve can be performed 
cheaply with a three term recurrence. The second consequence is that a LU factorization of He 
can be incrementally computed and this factorization without pivoting is known to be stable for 
positive definite matrices. The conjugate gradient method is a very elegant, sophisticated and 
powerful algorithm that exploits nicely all the above mentioned properties. It can be implemented 
through short recurrences that do not require to store the complete set of vector Ve leading to a 
very low memory consumption. Furthermore, the conjugate gradient enjoys a unique minimum 
norm property on the forward error that reads xe = sj:gtmn xexo+! ^/ Aro \ \\x — x*\\a where x* 
denotes the exact solution and || • \\a is the norm associated with A. In addition, it exists an 
upper bound on its convergence rate that reads (£ > 1) 

\\x e -x*\\ A <2. \\xo-x*\\ A . (6) 

We refer to [551 [Ml [HO] f° r an exhaustive and detailed exposure of CG and to [39] for a nice 
description of its history. 
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3 Non- Hermit ian matrices 

In this section we detail augmentation and deflation techniques in Krylov subspace methods 
when the coefficient matrix A is non-Hermitian. We specifically focus on minimum residual 
norm subspace methods and assume that a generalized Arnoldi relation © holds. We denote 
by so, tq = b — Axq the initial guess and residual vector respectively, and by Ve+i, He and Zg 
the matrices involved in this relation. With notation of Algorithm [1] r$ can be expressed as 
ra = Vt+i{c - H e y*). 

3.1 Augmented Krylov subspace methods 

We next discuss two possibilities to select the augmentation space and analyze the corresponding 
Krylov subspace methods. 

3.1.1 Augmentation with an arbitrary subspace 

Given a basis W = [w±, ■ ■ ■ ,Wk] of an augmentation subspace W, a slight modification in the 
Arnoldi procedure (Algorithm [2]) is used to deduce an orthogonal basis of Si defined in ((1} 
(see [2]). ^ consists of defining Zj (line 2 of Algorithm [5]) now as 

Zj = MjVj (1 < j < to) and Zj = Mj Wj- m (to < j < m + k). 

With this definition we finally obtain the generalized Arnoldi relation 

AZ m+ k — V m +k+iH m+ k 

where 

Z, m+k = [Mf V, M^v 2 , • • ■ , M^ +lWl ,M-\ 2 W2, • • ■ , M~l k+1 w k ], (7) 

Vm+k+l = [VI,V2,--- ,V m +k+l], (8) 

and H m+ k is a (to + k + 1) x (m + fc) upper Hessenberg matrix. Thus the residual minimization 
property is then deduced similarly as in FGMRES [S7]. Hence, the approximate solution from 
the affine space xq + range(Z m+ fc) can be written as 

Xm+k = + Z m+k y* 
with y* G ^(m+k) som tion of the residual norm minimization problem 

y* = argmin ||||r ||ei - H m+k y\\ 

j /e c< m+fc) 

(with ei designing here the first canonical vector of R( m+fc+1 )). In case of constant right pre- 
conditioning the main important property is that if any vector vjj is the solution of AM~ 1 Wj = 
Vi,l < i < to, then in general the exact solution of the original system © can be extracted 
from Sn\ see, e.g., [SSI Proposition 2.1]. We refer the reader to [2] for a discussion of possible 
choices for the augmented subspace W. Vectors obtained with either different iterative meth- 
ods or with different preconditioners can be incorporated in Z m quite easily. In block Krylov 
subspace methods we also mention that W consists of the sum of a few other Krylov subspaces 
generated with the same matrix but with different right-hand sides; see |14) for a discussion and 
numerical experiments on academic problems. A popular idea is to choose VV as an approximate 
invariant subspace associated with a specific part of the spectrum of A or AM -1 in case of fixed 
preconditioning. This is discussed next. 
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3.1.2 Augmentation with approximate invariant subspace 

A typical goal of augmentation is to add information about the problem into the search space 
that is only slowly revealed in the Krylov subspace itself. It is often known that eigenvalues of 
the (preconditioned) operator close to zero tend to slow down the convergence rate of the Krylov 
subspace methods |14j . Hence, augmentation based on approximate invariant subspaces made of 
eigenvectors corresponding to small in modulus eigenvalues of the (preconditioned) operator has 
been proposed; see, e.g.,[321 ESI [Ml [S3 an d references therein. 



Harmonic Ritz information In [55] Morgan has suggested to select W as an approximate 
invariant subspace and to update this subspace at the end of each cycle. Approximate spectral 
information is then required to define the augmentation space. This is usually obtained by 
computing harmonic Ritz pairs of A with respect to a certain subspace |141I56] . We present here 
a definition of a harmonic Ritz pair as given in [63, 75 . 

Definition 2 Harmonic Ritz pair. Consider a subspace U of C™. Given B G C nxn , 6 € C and 
y £ U, (6, y) is a harmonic Ritz pair of B with respect to IA if and only if 

By — 6y± BU 

or equivalently, for the canonical scalar product, 

Mw G Range(5 U) w H {By - 6y) = 0. 

We call the vector y a harmonic Ritz vector associated with the harmonic Ritz value 6. 

Based on the generalized Arnoldi relation ([3]), the augmentation procedure proposed in |37l 
Proposition 1] relies on the use of k harmonic Ritz vectors Y k = VeP k of AZeV e H with respect 
to Range (1^), where Y k € C nxk and P k = [pi, • • • ,p k ] G C exk . According to Definition!! the 
harmonic Ritz vector yj = Vepj then satisfies 

ZfA H (AZ tPj -e i Vip j ) - 0. (9) 

Using the generalized Arnoldi relation ([3]) we finally obtain the relation 

Hf H e Vj = eSfV&Vt Vj. (10) 

Since 



Hp = 



Hg G C 



where He G C £x£ is supposed to be nonsingular, the generalized eigenvalue problem is then 
equivalent to 

(He + h 2 e+he H- H eeeJ) yj = d ]Vr (11) 

This corresponds to a standard eigenvalue problem of dimension i only, where i is supposed 
to be much smaller than the problem dimension n. In consequence, the approximate spectral 
information based on Harmonic Ritz pair is quite inexpensive to compute. 
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GMRES augmented with approximate invariant subspace The augmentation 
space W based on approximate invariant information corresponding to Range(Yfc) is then used. 
The key point detailed next is to understand how to incorporate this information in a minimum 
residual norm subspace method such as GMRES. To do so, we recall a useful relation satisfied 
by the harmonic Ritz vectors P k £ C exk shown in |37[ Lemma 3.1] 



AZ e P k 



V l+ i 



Oixfe 

AZ e P k 



, c 



-Biy* 

[ViP k ,r 



diag(0i, ...,9k) 

diag(0i,...,0 fc ) 
aixfe 



(12) 
(13) 



where ro = Vg+i(c — Hgy*) and ai xk = [oti, ■ ■ ■ , oik] £ C lxfc . Next, the QR factorization of the 
(£ + 1) x (k + 1) matrix appearing on the right-hand side of relation (fT2]) is performed as 



Pk 

Olxfc 



, c - H e y* 



= QR 



(14) 



where Q £ c(^+i)x(fe+i) j ias orthonormal columns and R G Q(k+i)x(k+i) j s U pp er triangular, 
respectively. Then it can be shown that the relations 



AZ k 
V k H +1 V k+ i 
Range ([Y k ,r ]) 

hold with new matrices Z k , V k G C nxfc 



V k+ \H k , 
Ik+i, 

Range(Vfc+i), 
and H k G C< fc+1 ) xfc defined as 

Z k = Zi Qixk, 

V k+1 = Vt+x Q, 

Hk = Q H Hi Qtxk, 



(15) 
(16) 
(17) 

(18) 
(19) 
(20) 



where V^+i, Z( and Hi refer to matrices obtained at the end of the previous cycle; see [371 
Proposition 2]. With the augmentation subspace W = Range (Y k ), m Arnoldi steps with variable 
preconditioners and starting vector v k+ i are then carried out while maintaining orthogonality to 
V k leading to 



A [zfc+i, • • • , z m+k ] = [vk+i, ■ ■ ■ , v m+ k+i] H m and V m+k+1 V, 



m+k+l 



= 1 



m+k+l- 



We note that H m £ (£(m+i)xm - m U pp er Hessenberg. At the end of the new cycle this gives the 
generalized Arnoldi relation 



A [Zk, z k+ \, ■ ■ ■ , z m +k] — [V m +k+i] 





' H k ' 




Bkxm 






Omxfe 




H m 





i.e. 



AZ m+k — V m +k+iH m +k, 
H m+k £ C(" l + fe + 1 ) x (™+ fc ) and B kxr 



~>k x m 



results from the 



where V m + k +i £ C nx< - m+k+1 \ 
orthogonalization of [Az k +i, ■ • • , Az m + k +i] against V k +i- We note that H m+k is no more upper 
Hessenberg due to the leading dense {k + 1) x k submatrix H k . It is important to notice that 
the augmentation space varies at each restart since it is built from the search space available at 
the end of each previous cycle. The resulting algorithm can be viewed as an adaptive augmented 
Krylov subspace method. We refer the reader to [37J Sections 2 and 3] for the complete derivation 
of the method and additional comments on its computational cost. 
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Remarks and applications When the preconditioner is fixed, the previous algorithm pro- 
posed by Morgan [55] is known as GMRES with deflated restarting (GMRES-DR). Although the 
term "deflated" is used, we note that this algorithm does correspond to a GMRES method with 
an adaptive augmented basis without any explicit deflated matrix. The success of GMRES-DR 
has been demonstrated on many academic examples |54| and concrete applications such as in 
lattice QCD [T5I I29|. reservoir modeling (21115] or electromagnetism [37]. We refer the reader 
to (55J [55J for further comments on the algorithm and computational details. We note that 
GMRES with deflated restarting is equivalent to other augmented GMRES methods such as 
GMRES with eigenvectors [M] and implicitly restarted GMRES [55]. Most often the approxi- 
mate invariant subspace is chosen as the Harmonic Ritz pair corresponding to the smallest in 
modulus Harmonic Ritz values. Depending on the problem we note however that other specific 
part of the spectrum of the preconditioned operator can be targeted; see, e.g., |37l Section 4.2] 
for an application related to a wave propagation problem. 

3.2 Deflated Krylov subspace methods 

We next briefly describe minimal residual Krylov subspace methods based on deflation. We refer 
the reader to [501 I3T1 |4T)] for a recent excellent overview of deflated Krylov subspace methods in 
the Hermitian and non-Hermitian cases, where extensive bibliographical references and historical 
comments can be found. The general idea of deflation is to split the approximation space into 
two complementary subspaces such that the projected linear system, referred to as the deflated 
linear system, will be easier to solve iteratively than the original linear system ((2j). The fact that 
these subspaces can be chosen in different ways explains the huge literature on deflated Krylov 
subspace methods. The Krylov subspace method is then confined in one of this subspace, by 
projecting the initial residual into this space and by replacing A by its restriction to this space. 
If the projection operator is chosen properly the deflated linear system will be easier to solve 
iteratively than the original linear system (J5J . This property will be notably shown for Hermitian 
positive definite systems in Section 14.21 and can be extended to non-Hermitian situations with 
additional assumptions on A (see, e.g., [551 Section 2]). We first present a possible strategy 
based on orthogonal projection and then briefly discuss an extension based on oblique projection 
proposed in |40| . 

3.2.1 Deflation based on orthogonal projection 

We still denote by W a subspace of C™ of dimension fc, where k is assumed to be much smaller 
than the problem dimension n. We later denote by W <S C nxk a matrix whose columns form a 
basis of W so that W H A H AW is HPD (hence invertible). To simplify further developments, we 
introduce the matrices Qi,Pi,P 2 G C" x ™ defined respectively as 



We can easily show that Pi and P 2 are orthogonal projectors such that Pi projects onto (AW) 1 - 
along (AW), whereas P2 projects onto W' 1 along W. Furthermore we note that Pi is Hermitian 
and that AP2 = Pi A. The decomposition based on orthogonal projection reads as 



Qi 
Pi 

Pi 



AW(W H A H AW)- 1 W H A H , 

In — Ql, 

I n - W(W H A H AW)~ 1 W H A H A. 



(21) 
(22) 
(23) 



C 



n 



Hence, the solution x* of the original system @ can be written as 

x* = {I- P 2 )x* + P 2 x* = W(W H A H ' AWy 1 W H A H b + P 2 x*. 
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With this decomposition the original system ^ then simply becomes 

P x Ax = Pxb. (24) 

Although the deflated matrix Pi A is singular, the deflated linear system (j2~4")l is consistent so that 
it can be solved by an appropriate Krylov subspace method. Here we focus on the application of 
minimum residual Krylov subspace method based on GMRES to solve the deflated linear system 
(|2"4"|) . Hence, the search space of the Krylov subspace method applied to (f2"4"|) can be written as 

S m = IC m {PiA, Pir ), 

while the current approximation x m and the current residual f m = Pi (6 — Ax m ) at the end of 
the cycle satisfies the relations 

*£m ^ SCO ~l~ S m , 

PL(b-M m ) 1 PiAK^CPi^Piro). 

Since P±AW = 0, lX fe, Pi A is singular. Hence it is of paramount importance to analyze the 
possibilities of a breakdown when solving the deflated linear system In our context, when 

GMRES is used to solve the deflated linear system, this feature has been notably analyzed in [HH 
Section 3] based on theoretical results obtained by Brown and Walker [9] . We refer the reader to 
[4"01 Corollary 3] for conditions that characterize the possibility of breakdowns. It is worthwhile 
to note that a breakdown cannot occur if the condition 

Ker(P x A) n Range (Pi A) = {0} 

holds; see |301 Theorem 4.1]. This condition is notably satisfied if W is chosen as an exact A- 
invariant subspace, i.e., when AW = W since Ker[P\A) = W and Im(PiA) = W due to the 
nonsingularity of A. Once the solution of the deflated linear system is obtained, we deduce the 
approximation x m of the original system as 

x m = W{W H A H AW)- Y W H A H b + P 2 x m , 

and by construction we note that 

b - Ax m = Pi (6 - Ax m ), 

i.e., 

We refer to [55] for applications of deflated Krylov subspace methods with orthogonal pro- 
jection to linear systems with non-Hermitian matrices. As an illustration, a typical choice of 
subspaces is to choose the columns of W as right eigenvectors of A corresponding to eigenvalues 
of small absolute value. 

3.2.2 Deflation based on oblique projection 

We briefly mention a strategy based on oblique projection that is considered as more appropriate 
for the solution of non-Hermitian linear systems since the eigenspaces of A are in general not 
mutually orthogonal [JD] • As in Section 13.2.11 the search space Si will be decomposed into a 
direct sum of two subspaces. More precisely, the following decompositions into nonorthogonal 
complements are used 

C" = AW © VV^ = AW ® W^, 



Inria 



Deflation and augmentation techniques in Krylov lienar solvers 



13 



where W and W represent two subspaces of C" of dimension k respectively. As before, we denote 
by W 6 <C nxk (W E C nxk ) a matrix whose columns form a basis of W (W respectively). We 
assume that both matrices are chosen such that W H AW is nonsingular. The key idea is then to 
introduce the matrices Q 2 ,P 3 G C nx " defined as 

Q 2 = W(W H AW)~ 1 W H ', (25) 
P 3 = I n -W{W H AW^W" . (26) 

It is easy to show that Q2 and P 3 = I n — Q2 are projection operators; Q2 is the oblique projection 
onto (AW) along VV , while P 3 is the oblique projection onto W' L along (AW). Given these 
oblique projection operators, the deflated linear system is now defined as 

P 3 AP 3 x = P 3 b 

with fo = P 3 rQ. The use of a Krylov space solver is then now restricted to VV -1- . Hence, it can 
be shown that the deflated Krylov subspace method based on GMRES yields iterates x m at the 
end of the cycle of the form 

+ IC m (P 3 AP 3 ,P 3 r ) + W. 

This also implies the following relation for the residual |40j 

b-Ax m e r +AJC m {P 3 AP 3: P 3 r )+AW. 

We refer the reader to [301 Sections 5 and 6] for the mathematical aspects of deflated Krylov 
subspace methods based on oblique projection and to [201 Section 11] for an overview of partly 
related methods that only differ in the choice of the projection operators. A typical choice is to 
choose the columns of W as right eigenvectors of A and the columns of W as the corresponding 
left eigenvectors. We refer to [55] for an application of deflated Krylov subspace methods with 
oblique projection in the general non-Hermitian case. 

3.2.3 Deflation by preconditioning 

Finally, we note that deflation based on spectral approximate information can be used to con- 
struct nonsingular preconditioners that move small in modulus eigenvalues away from zero. Both 
Kharchenko and Yeremin J35] and Erhel et al. [23] have proposed GMRES algorithms with aug- 
mented basis and a nonsingular right preconditioner that move the small eigenvalues to a (mul- 
tiple) large eigenvalue. Baglama et al. [5] have proposed a left preconditioned GMRES method 
with similar effect. In [55] the main idea is to translate a group of small eigenvalues of A via 
low-rank projections of the form 

A = A(I n + Ul w?) ■ ■ ■ (I n + u fc wf ), 

where uj and Wj are the right and left eigenvectors associated with the eigenvalues to be trans- 
lated respectively. The restarted Krylov subspace method is now applied to the coefficient matrix 
A leading to an adaptive update of the preconditioner (performed at the end of each cycle) . We 
note that A can correspond to an already preconditioned operator, in such a case this strategy 
leads to a two-level preconditioning strategy that is found to be effective on real-life applications 
provided that the spectral information is computed accurately |12| . We also mention the exten- 
sion of this two-level preconditioning strategy to the case of sequences of linear systems (see, 
e.g., |36| where additional theoretical results and numerical experiments can be found). 
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3.3 Augmented and deflated Krylov subspace methods 

In the previous sections, we have described how either augmentation or deflation can be incorpo- 
rated into Krylov subspace methods of minimum residual norm type. We note that it is possible 
to combine simultaneously deflation and augmentation in a single Krylov subspace method. In 
such a setting, the search space of the Krylov subspace method is then decomposed as 

S e = W + K m (A,r ) 

where W is the augmentation space of dimension k, A refers to the deflated operator and fg to 
the deflated residual. As an illustration, we review the GCRO (Generalized Conjugate Residual 
with Orthogonalization) method due to de Sturler [16] . 



3.3.1 Equivalence between deflated and augmented methods 

In this section, we describe a general setting that helps us to understand the link between deflated 
and augmented minimal residual norm Krylov subspace methods. It has been first presented in 
|40| and we generalize this setting to the case of flexible methods. As discussed in Section 13.11 
the search space in augmented methods is of the form 

S e = W®K m (A,r ) 

where W is an augmentation subspace of dimension k. The approximation x m at the end of a 
given cycle can be written as 

x m =x Q + Z m y m + Ww m 

where y m e C m and w m G C fc . In the augmented Krylov subspace methods that we have 
considered, the residual r m satisfies a Petrov-Galerkin condition, i.e., r m J_ AS which leads to 
the two orthogonality conditions 

r m 1 AW and r m 1 AK m {A, r ). 

The first orthogonality condition r m _L AW leads to the relation 

(W H A H AW)w m = W H A H (r Q - AZ m y m ). 

To simplify notation we introduce the matrix Q3 S C" x ™ such that 

Q 3 = W{W H A H AW)~ 1 W H . 

We then deduce the following relations for the current approximation x m 

X m = {In - Q 3 A H A){X + ZmVm) + Q 3 A H b 7 (27) 

and for the current residual r rn 

r m = (In - AQ 3 A H ){r - AZ m y m ). (28) 
We then introduce the two matrices 

Pa = I n -Q 3 A H A, 
P 5 = I n -AQ 3 A H 
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where P4 G C lixn and P5 G C nxn . It is easy to show that both P4 and P5 are orthogonal 
projectors and that AP4 = P5 A If we define x m G C™ as x m = x + Z m y m then relations ([27)) 
and (|2"5)) become 

x m = P4X m + Q 3 A H b, (29) 
r m - P 5 (6-AS m ). (30) 

Finally the second orthogonality condition r m _L v4/C m (A,fo) can then be stated as 

r m = P 5 (b - Ax m ) _L AK m (P 5 A, P 5 r ). 

We summarize these developments in the following proposition (see [301 Theorem 2.2]). 

Proposition 1 The following two sets of conditions 

x m G x + W + lC m (A,r ), 

b - Ax m 1 (AW + AIC m {A,f )), 

and 

x m G Xo+K. m (A,r ), 

f m = P 5 (b - Ax m ) L AIC m {A,f ) 

are equivalent in the sense that 

x m = PiX m + Q?,A H b and r m = f m . (31) 

The first set of conditions corresponds to the standard augmentation approach described in 
Section 13.11 In this class of methods the augmentation space W is explicitly included in the 
search space S of the minimum residual Krylov subspace method and A = A, fo = ro. The 
second set of conditions corresponds to the standard deflation approach described in Section 
13.21 Indeed the iteration x m is first obtained such that the residual P$(b — Ax m ) satisfies the 
Petrov-Galerkin orthogonality condition. Then a correction is added such that r m = f m . Both 
approaches are found to be equivalent. They only differ in the way the augmentation subspace 
is treated (explicitly or implicitly). 

3.3.2 Methods based on augmentation and deflation 

Methods based on both augmentation and deflation have been introduced recently; see, e.g., 
[HI Uni HZl [ST] • We focus here on the Generalized Conjugate Residual with inner Orthogonalization 
(GCRO) [T5], which combines augmentation and deflation judiciously as detailed next. 

GCRO belongs to the family of inner-outer methods [21 Ch. 12] where the outer iteration is 
based on the Generalized Conjugate Residual method (GCR), a minimum residual norm Krylov 
subspace method proposed by Eisenstat, Elman and Schultz |22j while the inner part is based 
on GMRES respectively. Following the theoretical framework introduced in [21], GCR main- 
tains a correction subspace spanned by Range(Zfc) and an approximation subspace spanned by 
Range(Vfe), where Z^, Vk G C nxfc satisfy the relations 

AZ k = V k , 

v k H v k = I k . 
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The optimal solution of the minimization problem min \\b— Ax\\ over the subspace XQ + Kange(Z k ) 
is then found as x k = xq + Z k V k H tq . Consequently r k = b — Ax k satisfies 

r k = r - V k V h H r = (J„ - V k V h H )r , r k 1 Range(14). 

In [TB] de Sturler suggested that the inner iteration takes place in a subspace orthogonal to the 
outer Krylov subspace. In this inner iteration the following projected linear system is considered 

(In - V k V k H )Az = (J„ - V k V k H )r k = r fc . 

The inner iteration is then based on a deflated linear system with (J„ — V k V k ) as orthogonal 
projection. If a minimum residual norm subspace method is used in the inner iteration to 
solve this projected linear system approximately, the residual over both the inner and outer 
subspaces are minimized. Hence, augmentation is applied in the outer iteration and deflation 
in the inner part of the method. Numerical experiments (see, e.g., |16| and [271 Chapter 1]) 
indicate that the resulting method may perform better than other inner-outer methods (without 
orthogonalization) in some cases. 

We mention that the augmentation subspace can be based on spectral approximate invariant 
subspace information. This leads to the GCRO with deflated restarting method (GCRO-DR) 
|65| that uses Harmonic Ritz information to define the augmentation subspace as in Section 
13.11 This method has been further extended to accommodate variable preconditioning leading 
to the FGCRO-DR method [T3]. We also refer the reader to Q2] for additional comments 
on the computational cost of FGCRO-DR and a detailed comparison with the flexible variant 
of GMRES-DR. When a fixed right preconditioner is used, GMRES-DR and GCRO-DR are 
equivalent. When variable preconditioning is considered, it is however worthwhile to note that 
FGMRES-DR and FGCRO-DR are only equivalent if a certain collinearity condition given in 
[131 Theorem 3.6] is satisfied. 

In |17| de Sturler proposed to define an augmentation subspace based on information other 
than approximate spectral invariant subspace. At the end of each cycle, the strategy (named 
GCRO with optimal truncation (GCROT)) decides which part of the current global search 
subspace to keep to define the new augmentation subspace such that the smallest inner residual 
norm is obtained. This truncation is done by examining angles between subspaces and requires 
specification of six different parameters that affect the truncation. We refer to [17] for a complete 
derivation of the method and numerical experiments (see also [2T1 Section 4.5]). Finally we note 
that the extension of GCROT to the case of variable preconditioning has been proposed in [42] 
with application to aerodynamics. 

4 Hermitian positive definite matrices 

Similarly to unsymmetric problems both augmentation and deflation can be considered to speed- 
up the convergence of the conjugate gradient method, possibly in combination with precondi- 
tioning. However, contrarily to the previous methods based on Arnoldi basis construction, the 
conjugate gradient method relies on a short term recurrence and restarting mechanisms do not 
need to be implemented to control the memory consumption. Consequently the space used for 
augmentation or for deflation should be fully defined before starting the iteration for a given 
right-hand side. 

4.1 Augmented conjugate gradient methods 

As discussed in Section |3~T1 the search space in augmented methods Si = W © K m (A,ro) is a 
£ dimensional space (with i = m + k) where W is an augmentation subspace of dimension k 
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spanned by k linearly independent vectors W = \w\, ■ ■ ■ Wk\- In order to build a basis of this 
space a deflated Lanczos algorithm can be used, that consists in applying a standard Lanczos 
method starting from an unit vector v\ using the matrix 

B = A — AW(W H AW)- 1 W H A. 

Notice that W H AW is full rank since A is HPD. If v\ is orthogonal to W, deflated Lanczos 
builds a sequence of orthonormal vectors V m = [vi, ■ ■ ■ ,v m ] (VJ^AV m = I m ) that spans a space 
orthogonal to W, i.e., W H V m = Okxm- F° r the Krylov subpace part of St, xo is chosen so 
that ro = b — Axq _L W and v\ = ?*o/||ro||. That can be guaranteed by defining Xq from any 
X-i as xo = X-x + W(W H AW)~ 1 W H r-i. The augmented CG algorithm seeks for a solution 
xe — xq + WfJ,£ + V m y£ G Xq + W + K m {A, ro) with the Ritz-Galerkin condition rg = b — Axe _L 
(W + K m (A,ro)). Using the above described space and orthogonality condition, it is shown 
in [71], that the following properties (that are very similar and inherited from the classical CG) 
still hold. 

Proposition 2 The iterate Xj, the residual rj and the descent directions pj satisfy the following 
relations and properties 

• rj is collinear to Wj+i, that is, the residual vectors are orthogonal to each other, 

• the short term recurrences are satisfied: 

Xj = Xj-t + Ofj-lPj-l 

r = rj-! - a 3 -!Ap 3 -! 

P j = Tj+pj-iPj-i-WHj 

where ctj-! and f3j—! have the same expression as in classical CG and 

fij = (W H AW)- 1 W H Ar j , 



• the vectors pj are A-orthogonal to each other as well as A-orthogonal to all the Wj 's. 

Using theoretical results from [25], the following properties related to convergence rate and 
optimization property of the iterate are shown in [71] . 

Proposition 3 The approximate solution xg is such that 

• the convergence history exhibits an upper bound expression on the convergence rate similar 
to classical CG 



Ik(P* Xa AP w x a )-1 
\xt-x*\\ A <2- | V = = I \\xq-x*\\ a , (32) 



where «;(•) denotes the condition number of the matrix and P\^-l a is the A-orthogonal 
projection on W J ~ A . This projection is defined by P W ± A = I n — W(W H AW)~ 1 W H A. 

similarly to classical CG, the iterate complies with a minimum A-norm error on the search 
space xi = argmin || x — x* \\ a- 

xEx +W+K m (A,r ) 
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Augmenting using an invariant subspace Let (Ai = A m j n , • • • , X s = A max ) denote the s 
distinct eigenvalues of A ordered by increasing magnitude (i.e., values as they are real positive). 
The invariant subspace spanned by the k extreme (either largest or smallest) eigenvalues can 

be used in place of W to build the augmented space. Equation (j3"2")) shows that n(A) ~ max 

^min 

(that would appear in this bound for classical CG) is replaced either by max if the left most 

Afc+1 

X s — k 

part of the spectrum is used or by — if the right most part is used. Consequently if A m ; n <C 

Amin 

Xk+i (A s _fc <C A max ) the convergence of augmented CG should be significantly faster than the 
convergence of CG on the original system. 



4.2 Deflated Krylov subspace methods 

We next briefly describe CG variants based on deflation. As mentioned in Section l3~2l the general 
idea of deflation is to split the approximation space into two complementary subspaces. Similarly 
to the notation in the previous sections, we denote by W a subspace of C" of dimension k, where 
k is assumed to be much smaller than the problem dimension n. We later denote by W € C nxk 
a matrix whose columns form a basis of W. Because A is Hermitian positive definite, W H AW 
is also HPD and hence invertible. We can then define the following projector 

P 6 = I -W(W H AW)~ 1 W H A (33) 

that is an oblique projector along W {Pq is equal to P W ± A ). As in the non- Hermitian case, 
we decompose the solution x* = (I — Pq)x* + Pqx* and compute each component separately. 
In particular, (I - P 6 )x* = W(W H AW)- X W H Ax* = W(W H AW)- 1 W H b essentially reduces 
to the solution of a small k x k system. For the calculation of the second component Pqx* , 
it can be observed that AP 6 = P 6 H A so that AP 6 x* = P^ Ax* = P^b. Even though the 
matrix P§ A is Hermitian semi-definite positive of rank n — k (its nullspace is W), CG can 
still be used because the deflated linear system P^ Ax* = P§b is consistent jJS]. Furthermore, 
because the null space never enters the iteration, the corresponding zero eigenvalues do not 
influence the convergence |45j and we can define the effective condition number of the positive 
semidefinite matrix P^ 1 A, denoted K e ff(P^ A), as the ratio of its largest to smallest strictly 
positive eigenvalues. 

Once the linear system Pg Ax = P^b is solved, one just needs to apply Pq to this solution 
to compute the second component of the solution. This technique still requires the solution of a 
linear system of size n using the CG method, but is expected to be more effective if K e f /(P@ A) <g; 
k(A). We refer the reader to [28] for a discussion on the choice of W. 



Deflating using an invariant subspace If W defines an invariant subspace of A associated 
with extreme eigenvalues, the situation becomes much clearer. 

Let assume that W defines an invariant subspace associated with the smallest eigenvalues 
(Ai, Afe) of A. We have Pg AW = mX k so that P§ A has k zero eigenvalues. Because A is 
HPD, Z = W 1 -, the orthogonal complement of W (i.e., W H Z = so that P<f Z = Z) defines 
an invariant subspace associated with the eigenvalues Afc+i,...,A„ = A max - Therefore, we have 
AZ = ZB for some nonsingular B. Consequently we have P§ AZ = P^ 1 ZB = ZB so that Z is 
an invariant subspace of P§ A associated with the same eigenvalues Afc+i, A max . This shows 
that 

Keff (P 6 H A) = X 



A 



k+l 
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that indicates that deflating using an invariant subspace cancels the corresponding eigenvalues, 
leaving the rest of the spectrum unchanged. If Afc + i ^> Ai = A min the convergence of CG is 
significantly speeded-up. 

4.3 Deflation via preconditioning 

Using spectral information, it is possible to design preconditioners that enable to exhibit a 
condition number for the preconditioned matrix similar to K e ff{P^ A). 

Let W = [wi, . . . Wk] S C nxfc be the normalized eigenvectors of A associated with {Ai}t=i,„. fc 
the set of smallest eigenvalues. Let v be a real positive value. We can then define the precondi- 
tioner 

M def = i n + w(v(w H Awy 1 - h)W H . 

This preconditioner is such that M de f AW = vW and M de f Aw = Aw if W H w = (in particular 
any eigenvectors of A not in W), which shows that M de f moves the eigenvalues {At}t=i,. * to 
v and leaves the rest of the spectrum unchanged. If v = Xk+i, the condition number of the 
preconditioned matrix is the same as the one of the deflated matrix in the previous section. 

Furthermore we can define additive coarse space correction preconditioners inspired from 
domain decomposition techniques. They lead to preconditioned matrices with similar condition 
number as well. We then define 

M coarse = J n + V W {W H AW^W" . 

This preconditioner is such that M coarse Awi = {v + X^Wi and M coarse Aw = Aw if W H w = 0. 
That is, the eigenvalues {Ai}t=i ... k are shifted to v + Ai, while the rest of the spectrum is 
unchanged. If it exists v so that Afc + i < A m ; n + v < Xk + v < A max , the preconditioned matrix 
would have again the same condition number as the one of the deflated system K e ff- 

We refer the reader to [3S] for an analysis of the condition number of this class of precondi- 
tioners when approximated spectral information is used. We also refer to |78| and the references 
therein for the exposure of various preconditioning techniques that can be defined using various 
combinations of these building box components. 

5 Linear systems with multiple right-hand sides given in 
sequence 

Although our primary focus is the solution of a single linear system with preconditioned Krylov 
subspace methods, it is however possible to include deflation and augmentation in a broader 
setting. Indeed in many applications in computational science and engineering, linear systems 
with multiple right-hand sides have to be solved. More precisely we are interested in solving 
a sequence of linear systems defined as A 1 x l = b l where both the non-Hermitian matrix A 1 £ 
C nxn and the right-hand side b l e C" may change from one system to the next, and the linear 
systems may typically not be available simultaneously. If we consider a sequence of identical 
or slowly changing matrices, Krylov subspace methods based on augmentation and deflation 
are appropriate since subspace recycling is then possible. The key idea is to extract relevant 
information (e.g. approximate invariant subspace but not only) while solving a given system, 
and then to use this information to further accelerate the convergence of the Krylov subspace 
method for the subsequent linear systems. At this point, augmented or deflated Krylov subspace 
methods of Sections 13.11 13.21 and 13.31 can then be used. We refer the reader to |64l Chapter 
3] for a detailed analysis of subspace recycling in the non-Hermitian case and to |65| where 
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the GCRO method augmented with approximate spectral information is shown to be efficient 
on applications related to fatigue and fracture of engineering components, electronic structure 
calculation and quantum chromodynamics; see also |47| for an application in optical tomography. 
Recent applications are related to model reduction [7] (see also [T] for recycling methods based 
on BiCG). 

For HPD matrices, if a sequence of linear systems with the same matrix but different right- 
hand sides has to be solved different alternatives can be considered to define the space to augment 
the search space from one solve to the next. In |25j . an approach based on harmonic Ritz values is 
described that might be implemented using only the first m > k steps of augmented CG iteration. 
Still to reduce the memory footprint of the eigenvector calculation, in [76] a thick-restart Lanczos 
is embedded in the CG iterations to extract accurate spectral information. 

6 Conclusions and prospectives 

We have briefly reviewed the main features and mathematical properties of augmented and 
deflated Krylov subspace methods for the solution of certain linear systems of equations where 
the coefficient matrix was either non-Hermitian or Hermitian positive definite. These increasingly 
popular procedures combined with preconditioning have been shown effective on a wide range 
of applications in computational science and engineering as mentioned in this paper. We are 
certainly aware that this brief overview is far from being complete. Results related to two-sided 
Krylov subspace methods in the non-Hermitian case or the treatment of the Hermitian indefinite 
case are indeed missing; see, e.g., [1, 30 ( 3f , 83J for additional comments and references. Similarly, 
the solution of linear systems with multiple right-hand sides given at once has not been covered. 
For such a class of problems, augmented and deflated block Krylov subspace methods have been 
studied (see, e.g., [571 I84j ) and their efficiency has been proved on realistic applications. Finally 
we would like to mention that algebraic connections between deflation, multigrid and domain 
decomposition have been made in recent papers 051 E3 EFJ ■ 

Concerning implementation aspects, some augmentation and deflation procedures are already 
present in the main software projects such as either PETS(0orTrilino<J3 for the solution of large- 
scale, complex multi-physics engineering and scientific problems. More precisely, in its scalable 
linear equation solvers (KSP) component, PETSc includes an algorithm described in [53], while 
the Belos package in Trilinos notably proposes an augmented and deflated approach based on 
GCRO-DR [65]. Most likely there will be a growing effort to incorporate augmented and deflated 
Krylov subspace methods in such libraries in a near future. Finally designing variants or new 
Krylov subspace methods for the next generation of massively parallel computing platforms is 
currently a topic of active research in the numerical linear algebra community; see [3311341 H5] for 
algorithms, comments and references. Thus in a near future it is highly probable that variants 
of augmented and deflated Krylov subspace methods will be proposed as well. 
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